rm(list = ls())
library(ggplot2)
library(matlab)
alpha <- 0.5
beta <- 0.25
gm <- 0.25
dlt <- 0.5
theta <- 1
eta <- 1

# initial value
z1 <- 1
z2 <- 1

# steady state for initial value
A <- matrix(c(-alpha,beta, gm,-dlt),2,byrow = T)
B <- matrix(c(theta,0,0,eta),2,byrow = T)
Z1 <- matrix(c(z1,z2),2)
barX1 <- -(solve(A) %*% B %*% Z1)

# final value and steady state
z1 <- 2
z2 <- 1
Z2 <- matrix(c(z1,z2),2) # 注意z的大小写的含义是不同的
barX2 <- -(solve(A) %*% B %*% Z2)

# transition dynamics:从一个稳态到另一个稳态，中间的转移动态，也叫脉冲响应分析
xt <- matrix(NA,30,2)
xt[1,] <- t(barX1[,1])
for (i in 2:nrow(xt)) {
  xt[i,] <- t((A + eye(2)) %*% matrix(xt[i-1,],2) + B %*% Z2)
}

# draw
ggplot(data = as.data.frame(xt), aes(x = 1:nrow(xt), y = V1)) + geom_line() +
  geom_line(aes(y = V2), linetype = 2) + labs(x = 't', y = 'x') +
  theme_bw()
# ggsave('../arm_irf.pdf')


# change parameters: sensitivity analysis: 全局稳定
alpha <- 0.7
A <- matrix(c(-alpha,beta, gm,-dlt),2,byrow = T)
eigen(A)$values
-(solve(A) %*% B %*% matrix(c(1,1),2))

xt <- matrix(NA,30,2)
xt[1,] <- t(barX1[,1])
for (i in 2:nrow(xt)) {
  xt[i,] <- t((A + eye(2)) %*% matrix(xt[i-1,],2) + B %*% Z1)
}

# draw
ggplot(data = as.data.frame(xt), aes(x = 1:nrow(xt), y = V1)) + geom_line() +
  geom_line(aes(y = V2), linetype = 2) + labs(x = 't', y = 'x') +
  theme_bw()
# ggsave('../arm_sen1.pdf')

# change parameters: sensitivity analysis: 鞍点稳定
alpha <- 0.25
beta <- 0.5
gm <- 0.5
dlt <- 0.25

A <- matrix(c(-alpha,beta, gm,-dlt),2,byrow = T)
lamd <- eigen(A)$values
-(solve(A) %*% B %*% matrix(c(-1,-1),2))

# 鞍点稳定下，脉冲响应冲击:z1 = -0.5
barX3 <- -(solve(A) %*% B %*% matrix(c(-0.5,-1),2)) # steady status

xt <- matrix(NA,30,2)
xt[1,] <- t(barX1[,1])
for (i in 2:nrow(xt)) {
  if (i == 2){
    # 调整到鞍点路径
    xt[i,] <- t((A + eye(2)) %*% matrix(xt[i-1,],2) + B %*% matrix(c(-1,-1),2))
    xt[i,1] <- beta/(alpha+min(lamd))*xt[i,2] + theta/(alpha+min(lamd))*(-0.5)+
      min(lamd)/(alpha + min(lamd))*barX3[1]
  }else {
    # 根据鞍点路径运动
    xt[i,] <- t(min(lamd) * (matrix(xt[i-1,],2)-barX3) + matrix(xt[i-1,],2))
  }
}

# draw
ggplot(data = as.data.frame(xt), aes(x = 1:nrow(xt), y = V1)) + geom_line() +
  geom_line(aes(y = V2), linetype = 2) + labs(x = 't', y = 'x') +
  theme_bw()
# ggsave('../arm_irf2.pdf')
